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METHOD  FOR  DESIGN  ROTATION 


1.  INTRODUCTION 

Response  surface  designs  are  used  to  develop  a  polynomial  model  relating  a 
dependent  variable  y  to  a  set  of  quantitative  independent  variables  Zj,  Zg, . . . , zfc. 
Each  row  of  the  Nxk  design  matrix  X  is  a  point  in  k- dimensional  space.  The 
dependent  variable  y  is  measured  at  each  of  the  design  points  and  the  coefficients  of 
the  polynomial  model  are  estimated  from  the  data  by  the  least-squares  formula 

b  =  (M'M)-1M'y,  (1) 

where  b  is  the  estimated  coefficient  vector  and  M  is  the  Nxp  expanded  design 
matrix,  or  model  matrix,  that  has  one  column  for  each  coefficient  of  the  polynomial 
model. 

Many  designs  have  been  proposed  for  the  estimation  of  the  coefficients  of 
first-  and  second- degree  polynomial  models.  The  designs  are  centered  at  the  origin 
and  scaled  so  that  the  design  levels  are  convenient  numbers — such  as  integers, 
simple  fractions,  and  their  square  roots;  a  change  of  origin  and  scale  makes  the 
designs  applicable  to  applied  problems.  Because  a  response  surface  design  is  a  set  of 
points  in  i- dimensional  space,  a  rotation  of  the  coordinate  axes  will  produce  a  new 
design  matrix  that  represents  the  same  geometric  figure  as  the  original  design.  An 
elementary  result  in  linear  algebra  is  that  multiplication  by  an  orthogonal  matrix  P 
corresponds  to  a  rotation  of  the  coordinate  axes  if  the  determinant  of  P  is  1,  and  to 
a  rotation  followed  by  a  reflection — a  change  of  sign  or  direction  of  one  axis — if  the 
determinant  of  P  is  —1.  Hence  there  is  no  theoretical  problem  in  obtaining 
rotations  of  a  design;  there  are  only  the  practical  questions  Why  rotate  a  design? 
And  how  can  a  useful  orthogonal  matrix  be  selected  to  rotate  the  design? 


2.  MOTIVATION 

The  well  known  D-  and  G-efficiencies  of  first-  and  second-order  response 
surface  designs  are  not  affected  by  a  rotation  of  the  coordinate  axes — see  the  proof 
in  the  Appendix.  One  must  look  at  other  criteria  to  see  any  benefit  from  design 
rotation.  Two  criteria  examined  here  are  the  size  of  the  design  region  and  the 
symmetry  of  the  design  matrix.  The  starting  point  for  the  theory  of  response 
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surface  designs  is  the  specification  of  the  design  region;  it  is  usually  taken  to  be  a 
cube  or  sphere  of  a  specified  size  (in  coded  units).  For  spherical  designs,  the  ranges 
of  the  coded  design  factors  are  not  considered  important,  and  designs  are  compared 
by  scaling  them  to  the  same  diameter.  In  this  report,  I  take  a  different  approach  to 
developing  a  response  surface  design. 

Suppose  an  experimenter  wishes  to  vary  cooking  time  from  40  to  60  minutes, 
and  cooking  temperature  from  180  to  200  degrees.  I  assume  that,  regardless  of  the 
design  used,  the  range  of  the  coded  design  factors  must  be  scaled  to  the  ranges  of 
the  experimental  factors  so  that  cooking  time  varies  from  40  to  60  minutes,  and 
temperature  varies  from  180  to  200  degrees.  Such  a  scaling  of  the  coded  design 
factors  to  the  experimental  factors  is  necessary  to  give  the  experimenter  a  design 
with  the  specified  ranges. 

The  specification  of  ranges  for  the  experimental  factors  does  not  imply  a 
cuboidal  region.  The  midpoints  of  the  ranges  of  the  experimental  factors  usually 
represent  the  current  process  settings  or  a  guess  of  the  factor  settings  that  will 
optimize  y.  The  exploration  of  the  region  around  the  center  point  of  the  design 
should  be  symmetrical  (that  is,  a  spherical  design)  unless  the  experimenter  specifies 
that  some  directions  are  more  important  than  others.  A  cuboidal  region  would  be 
appropriate,  for  example,  if  the  experimenter  specifies  that  the  directions  toward 
the  comers  of  the  cube  are  more  important  than  the  directions  toward  the  centers 
of  the  faces  of  the  cube. 

Two  examples  will  show  how  the  size  of  the  design  region  and  symmetry  of 
design  matrix  can  be  changed  by  a  rotation  of  the  design.  The  first  example 
involves  a  rotation  and  scaling  of  the  central  composite  design  for  two  factors.  The 
central  composite  designs  were  developed  by  Box  and  Wilson  (1951)  and  consist  of 
(a)  center  points  at  the  origin,  (b)  axial  or  star  points,  which  can  be  represented  in 
a  shorthand  notation  by  the  permutations  of  (±1, 0,  •  •  •  ,0),  and  (c)  factorial 
points,  which  are  a  two -level  factorial  or  fractional  factorial  design  of  at  least 
resolution  V  and  can  be  represented  by  (±a-1,  ±a_1,  ,±a-1).  Typical 

choices  for  a  are  the  square  root  of  fc,  which  places  the  noncentral  design  points  on 
the  surface  of  a  sphere,  and  the  fourth  root  of  the  number  of  factorial  points,  which 
makes  the  contours  of  constant  prediction  variance  spherical.  Figure  1  shows  a 
rotation  and  scaling  of  the  two-factor  central  composite  design.  The  effect  of  the 
rotation  and  scaling  is  to  make  the  design  cover  a  larger  region  (the  large  outer 
circle  in  Fig.  l)  when  the  rotations  are  scaled  to  cover  the  same  range.  Table  1 
gives  the  coded  design  settings  of  the  central  composite  design,  the  rotated  design, 
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Temperature 


Time 


Figure  1.  A  Rotation  and  Scaling  of  the  Two-Factor  Central  Composite  Design. 
The  design  points  before  rotation  are  indicated  by  small  squares.  The  design  is  fust 
rotated  (solid  dots),  then  scaled  to  cover  the  same  range  as  the  design  before 
rotation  (hollow  dots). 


and  the  rotated  and  scaled  design.  The  advantage  of  the  larger  region  of  the 
rotated  and  scaled  design  is  that  it  gives  a  better  estimates  of  the  model  coefficients 
when  the  designs  are  applied  in  the  usual  manner.  Table  2  shows  that  the  larger 
experimental  region  of  the  rotated  and  scaled  design  in  Table  1  yields  better 
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Table  1.  Two  Factor  Central  Composite  Design 


Coordinates 

'...Rotated 

...and  Scaled 

Point 

*i 

Block 

1 

-.707  -.707 

-.924 

-.383 

-1.  -.414 

1 

2 

.707  -.707 

.383 

-.924 

.414  -1. 

1 

3 

-.707  .707 

-.383 

.924 

-.414  1. 

1 

4 

.707  .707 

.924 

.383 

1.  .414 

1 

5 

0.  0. 

0. 

0. 

0.  0. 

1 

6 

-1.  0. 

-.924 

.383 

-1.  .414 

2 

7 

1.  0. 

.924 

-.383 

1.  -.414 

2 

8 

0.  -1. 

-.383 

-.924 

-.414  -1. 

2 

9 

0.  1. 

.383 

.924 

.414  1. 

2 

10 

0.  0. 

0. 

0. 

0.  0. 

2 

Table  2.  Variance  of  Estimated  Coefficients 

Standard 

Rotated 

Model  Term 

Orientation 

and  Scaled 

constant 

.5 

<ra 

.5  a2 

linear 

.28 

<ra 

.213  a2 

quadratic 

.875  a 2 

.637  a2 

interaction 

1.0 

<ra 

.729  a2 

estimates  of  the  model  coefficients.  Scaling  the  design  to  cover  the  larger  region 
increases  the  determinant  of  M'M  from  256  to  909. 

The  second  example  is  a  rotation  and  scaling  of  the  hexagon  design.  Table  3 
gives  the  coordinates  of  the  design  both  before  and  after  the  rotation  and  scaling; 
Figure  2  shows  the  rotation  and  scaling. 


Table  3.  Hexagon  Design 
Coordinates  Rotated  and  Scaled 


*2 


Figure  2.  Two  Rotations  and  Scalings  of  the  Hexagon  Design.  The  hexagon  design 
in  its  usual  orientation  is  indicated  by  small  squares.  After  rotation  by  7r/12 
radians  and  scaling  the  factor  levels  to  the  range  —1  to  1,  the  design  points  are 
indicated  by  a  solid  dots  for  a  clockwise  rotation  of  the  design  and  by  hollow  dots 
for  a  counterclockwise  rotation  of  the  design.  Notice  that  the  clockwise  and 
counterclockwise  rotations  are  reflections  across  the  X-axis;  multiplying  the  second 
coordinate  by  —1  will  change  one  rotation  into  the  other. 

In  the  standard  orientation  of  the  hexagon  design,  the  first  factor  covers  the 
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interval  (—1,1)  and  the  second  factor  covers  the  interval  (—.866,. 866).  In  the  new 
orientations,  both  factors  cover  the  interval  (—1,1);  the  rotated  designs  also  cover  a 
slightly  larger  region  than  the  original  design. 

The  unequal  ranges  of  the  factors  of  the  hexagon  design  in  its  usual 
orientation  create  a  problem  in  applying  the  design.  If  the  levels  —.866  and  .866  of 
the  second  factor  are  scaled  to  180  and  200  degrees  for  the  two-factor  cooking 
example,  then  the  design  is  no  longer  a  spherical  design,  but  an  ellipsoidal  design 
based  on  an  irregular  hexagon.  This  method  of  applying  the  design  essentially 
ignores  the  levels  —.866  and  .866  of  the  second  factor;  they  could  be  replaced  by  -1 
and  1  without  affecting  the  levels  of  cooking  temperature  used  by  the  experimenter. 

An  alternative  method  of  applying  the  design  is  to  associate  the  diameter  of 
the  circle  with  the  ranges  of  the  experimental  factors.  The  design  diameter  of  two 
units  would  therefore  correspond  to  the  full  range  of  cooking  temperature  (180  to 
200  degrees).  Because  the  second  design  factor  varies  from  —.866  to  .866,  cooking 
temperature  would  be  varied  from  181.3  degrees  to  198.7  degrees  in  the  experiment. 
This  latter  method  of  applying  the  design  is  consistent  with  the  theory  of  optimal 
designs,  but  tends  to  produce  counter-productive  consultant-client  discussions. 

The  new  orientations  of  the  hexagon  design  give  both  factors  the  range  —1  to 
1,  so  they  do  not  create  a  problem  in  associating  the  range  of  the  coded  design 
factors  to  the  ranges  of  the  experimental  factors.  The  price  paid  for  the  symmetric 
orientations  of  the  hexagon  design  is  that  the  design  now  has  seven  levels  of  each 
factor.  An  alternative  is  to  use  the  central  composite  design,  which  has  only  five 
levels,  but  requires  more  experimental  runs. 

These  examples  suggest  that  both  the  practical  efficiency  and  the  correctness 
of  application  (or  the  client-friendliness)  of  a  response  surface  design  may  be 
affected  by  a  rotation  of  the  design.  So  how  can  desirable  rotations  be  obtained? 
There  is  little  information  in  the  literature  on  design  rotation.  Box  and  Behnken 
(1960)  note  that  their  four-factor  design  is  a  rotation  of  the  four-factor  central 
composite  design  and  give  the  orthogonal  matrix  that  yields  the  rotation,  but  they 
do  not  discuss  how  the  orthogonal  matrix  was  found. 

Doehlert  and  Klee  (1972)  gave  rotations  of  Doehlert’s  (1970)  uniform  shell 
designs  that  minimize  the  number  of  levels  of  the  factors.  A  technique  used  by 
Doehlert  and  Klee  (1972)  was  to  start  with  a  known  orthogonal  matrix  of  simple 
form  and  then  augment  the  matrix  with  additional  rows  and  columns  to  get  a 
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matrix  that  yields  the  minimum  number  of  factor  levels  for  the  next  higher 
dimensional  design.  The  augmentation  process  can  be  applied  to  other  sequences  of 
designs,  but  a  useful  orthogonal  matrix  is  required  to  start  the  process.  Crosier 
(1991)  gave  symmetric  orientations  of  response  surface  designs  for  mixture 
experiments.  No  orthogonal  matrices  were  given,  although  there  must  be  matrices 
that  map,  for  example,  the  three-dimensional  central  composite  design  to  the 
mixture  space  for  four  components. 

In  summary,  the  existing  literature  provides  little  guidance  for  rotating 
response  surface  designs.  The  only  helpful  idea  is  to  augment  an  orthogonal  matrix 
to  serve  the  next  higher  dimension  in  a  sequence  of  designs  that  have  the  same 
construction.  The  development  of  useful  orthogonal  matrices  to  start  the  process  is 
an  intuitive  process.  The  numerical  techniques  developed  in  this  report  are  offered 
as  a  supplement  to  the  intuitive  approach  to  obtaining  design  rotations.  Neither 
approach  is  guaranteed  to  find  the  best  solution.  Doehlert  and  Klee  (1972),  for 
example,  missed  the  five- level  orientation  of  the  eight-factor  uniform  shell  design 
given  in  Table  4.  The  design  in  Table  4  can  be  orthogonally  blocked:  the  points  in 
Table  4  are  one  block  and  their  negatives  are  the  other  block;  the  same  number  of 
center  points  must  be  used  with  each  block.  The  design  in  Table  4  is  one  design  in 
the  class  of  uniform  shell  designs;  the  symmetric  orientation  of  the  uniform  shell 
designs  was  found  using  the  numerical  techniques  described  in  the  next  section. 


3.  AN  APPROACH  TO  DESIGN  ROTATION 


3.1  Parameterization  of  Orthogonal  Matrices 

The  rotation  of  a  response  surface  design  can  be  written  in  matrix  notation  as 

W  =  XG  (2) 

where  W  is  the  rotated  design  matrix,  X  is  the  Nxk  design  matrix,  and  G  is  a  kxk 
orthogonal  matrix  whose  determinant  is  one. 

To  optimize  some  criterion  for  W  as  a  function  of  G,  it  is  necessary  to 
parameterize  G  as  a  vector  0  whose  elements  can  be  varied  independently  of  each 
other.  The  k2  elements  of  G  cannot  be  varied  independently  because  G  must 
remain  an  orthogonal  matrix.  G  can  be  expressed  as  the  product  of  k(k— 1)/2 
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Table  4.  Rotated  Eight -Factor  Uniform  Shell  Design 


Point 

*» 

** 

*» 

1 

-1. 

1. 

0. 

0. 

0. 

0. 

0. 

0. 

2 

0. 

-1. 

1. 

0. 

0. 

0. 

0. 

0. 

3 

0. 

0. 

-1. 

1. 

0. 

0. 

0. 

0. 

4 

0. 

0. 

0. 

-1. 

1. 

0. 

0. 

0. 

5 

0. 

0. 

0. 

0. 

-1. 

1. 

0. 

0. 

6 

0. 

0. 

0. 

0. 

0. 

-1. 

1. 

0. 

7 

0. 

0. 

0. 

0. 

0. 

0. 

-1. 

1. 

8 

1. 

0. 

0. 

0. 

0. 

0. 

0. 

-1. 

9 

-1. 

0. 

1. 

0. 

0. 

0. 

0. 

0. 

10 

0. 

-1. 

0. 

1. 

0. 

0. 

0. 

0. 

11 

0. 

0. 

-1. 

0. 

1. 

0. 

0. 

0. 

12 

0. 

0. 

0. 

-1. 

0. 

1. 

0. 

0. 

13 

0. 

0. 

0. 

0. 

-1. 

0. 

1. 

0. 

14 

0. 

0. 

0. 

0. 

0. 

-1. 

0. 

1. 

15 

1. 

0. 

0. 

0. 

0. 

0. 

-1. 

0. 

16 

0. 

1. 

0. 

0. 

0. 

0. 

0. 

-1. 

17 

-1. 

0. 

0. 

1. 

0. 

0. 

0. 

0. 

18 

0. 

-1. 

0. 

0. 

1. 

0. 

0. 

0. 

19 

0. 

0. 

-1. 

0. 

0. 

1. 

0. 

0. 

20 

0. 

0. 

0. 

-1. 

0. 

0. 

1. 

0. 

21 

0. 

0. 

0. 

0. 

-1. 

0. 

0. 

1. 

22 

1. 

0. 

0. 

0. 

0. 

-1. 

0. 

0. 

23 

0. 

1. 

0. 

0. 

0. 

0. 

-1. 

0. 

24 

0. 

0. 

1. 

0. 

0. 

0. 

0. 

-1. 

25 

-1. 

0. 

0. 

0. 

1. 

0. 

0. 

0. 

26 

0. 

-1. 

0. 

0. 

0. 

1. 

0. 

0. 

27 

0. 

0. 

-1. 

0. 

0. 

0. 

1. 

0. 

28 

0. 

0. 

0. 

-1. 

0. 

0. 

0. 

1. 

29 

.5 

-.5 

-.5 

-.5 

-.5 

-.5 

-.5 

-.5 

30 

-.5 

.5 

-.5 

-.5 

-.5 

-.5 

-.5 

-.5 

31 

-.5 

-.5 

.5 

-.5 

-.5 

-.5 

-.5 

-.5 

32 

-.5 

-.5 

-.5 

.5 

-.5 

-.5 

-.5 

-.5 

33 

.5 

.5 

.5 

.5 

-.5 

.5 

.5 

.5 

34 

.5 

.5 

.5 

.5 

.5 

-.5 

.5 

.5 

35 

.5 

.5 

.5 

.5 

.5 

.5 

-.5 

.5 

36 

.5 

.5 

.5 

.5 

.5 

.5 

.5 

-.5 

NOTE:  Add  negatives  and  center  points. 


matrices,  each  of  which  represents  a  rotation  in  a  plane  defined  by  two  coordinate 
axes.  All  pairs  of  axes  are  included  in  this  factorization  of  G,  which  can  be  written 

G  =  (G12Gi3  •  •  •  G1i!)(G23  •  •  •  Gj*)  •  •  •  (G*_ift),  (3) 

where  the  matrix  G,y  is  formed  by  replacing  four  elements  of  a  kxk  identity  matrix; 
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the  »,»  and  jj  elements  are  cos  (0tJj ,  the  ij  element  is  —sin  (fly)  and  the  j,i  element 
is  8tn(0ij).  The  6,y  matrices  represent  a  clockwise  rotation  of  the  design — that  is, 
the  design  is  rotated  from  the  positive  x y  axis  toward  the  positive  z,  axis.  The 
representation  of  the  matrix  6  as  the  product  of  k(k~l)/2  matrices  is  a  useful 
conceptual  tool,  but  it  is  not  necessary  to  actually  form  and  multiply  all  the  G,;'s  to 
obtain  G.  G  can  be  obtained  from  an  identity  matrix  in  k(k-l)/2  stages  that 
require  only  a  linear  transform  of  two  columns;  each  stage  corresponds  to  post¬ 
multiplication  by  a  G|;  matrix.  Post -multiplication  by  G,;  transforms  the  ith  and 
jth  columns  by  the  equations 

9  *nt  =  cos  (0,y)  gni  +  sin  (0t;)  gnj  (4) 

and 

9*nj  =  ~*in  (*»/)  9ni  +  cos  (Oij)  gnj  (5) 

where  n=l,2, . . .  ,k  is  the  row  index  and  g*n{  and  g*nj  are  the  values  in  the  ith  and 
jth  columns  after  the  transformation.  In  programming  equations  (4)  and  (5),  it  is 
necessary  to  put  the  value  of  g*n{  into  a  temporary  storage  location  so  that  gni  is 
still  available  when  equation  (5)  is  implemented;  it  is  also  customary  to  take  the 
calculation  of  cos{Bi})  and  sin  (fly  )  outside  the  loop  over  the  rows  of  G. 

Because  G  is  an  orthogonal  matrix  with  determinant  1,  equation  (2)  does  not 
allow  for  rotations  followed  by  a  reflection.  Sign  changes  of  the  factors  do  not  affect 
the  criteria  proposed  below  for  the  desirability  of  a  rotation.  Hence,  the  method 
developed  in  this  section  is  able  to  find  a  good  rotation  of  a  design,  although  it  may 
be  necessary  to  change  the  sign  of  a  factor  to  put  the  rotated  design  in  a  convenient 
form. 


3.2  Criteria  for  Selecting  a  Rotation 

The  examples  in  Section  2  show  what  properties  of  a  design  may  be  changed 
by  rotation  and  therefore  what  criteria  need  to  be  specified  for  the  objective 
function  of  the  multivariate  optimization  procedure.  Rotating  a  design  does  not 
change  its  diameter,  so  the  range  of  the  coded  factors  is  used  to  measure  the  size  of 
the  design.  As  can  be  seen  in  Table  1,  rotation  of  a  design  should  minimize  the 
range  of  the  coded  factors  so  that  rescaling  the  design  will  increase  the  size  of  the 
experimental  region.  Because  most  rotations  are  likely  to  have  different  ranges  for 
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the  factors,  the  largest  range  (#max)  will  be  taken  as  the  first  criterion  to  be 
minimized. 

The  second  example  in  Section  2  showed  that  equal  ranges  for  the  factors 
might  be  obtained  by  rotation.  To  obtain  this  goal,  the  difference  between  the 
largest  and  smallest  range 


Ail  —  Umax  Rmin  (®) 

may  be  minimized  by  the  optimization  procedure.  The  symmetric  orientation  of 
the  uniform  shell  designs  was  found  by  using  the  sum  of  the  two  criteria  R  maT  and 
AR  as  the  objective  function  for  the  multivariate  optimization.  Although  this 
success  showed  that  these  criteria  could  be  effective,  additional  criteria  were 
developed  to  indicate  other  properties  of  symmetric  designs. 

Design  factors  can  have  equal  ranges  without  being  defined  over  the  same 
interval;  for  example,  the  first  factor  may  cover  the  interval  —1  to  2,  whereas  the 
second  factor  covers  the  interval  —2  to  1.  It  is  also  possible  that  the  factors  may 
have  the  same  range,  but  the  range  is  asymmetric  about  zero,  the  center  of  the 
design.  For  example,  the  factors  may  all  cover  the  interval  -1  to  2,  even  though 
the  design  is  centered  at  the  origin.  To  develop  a  test  for  symmetric  ranges,  let  L; 
be  the  minimum  value  in  the  jth  column  of  the  rotated  design  matrix  W  and  let  Hj 
be  the  maximum  value  in  the  jth  column  of  W.  A  criterion  for  obtaining 
symmetric  ranges  is 


Rtyrn  —  E  I L)  Hi  I  •  (^) 

;=i 

This  criterion  will  be  zero  if  L}  =  —Fj  for  j  =  1,  2,  ...»  k — that  is,  if  all  the  factors 
have  symmetric  ranges. 

In  addition  to  having  the  same,  symmetric  range  for  all  factors,  a  symmetric 
design  should  have  the  same  set  of  factor  levels  for  all  factors.  An  exact  test  for  the 
same  set  of  factor  levels  for  all  factors  would  be  cumbersome,  so  the  following 
indicator  was  developed  to  search  for  designs  with  the  same  set  of  levels  for  every 
factor.  Find  the  sum  of  absolute  values  (SAV)  of  the  factor  levels  for  each  factor; 
take  the  difference  between  the  largest  and  smallest  of  these  sums.  The  difference 
will  be  zero  if  the  factors  have  the  same  set  of  levels.  The  quantity  to  be  minimized 
by  the  optimization  procedure  is 
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(8) 


N  N 

ASA  V  =  max£  |  ti/y  |  -  minj]  |  u>y  | . 

i  t=i  J  i=i 

The  criteria  developed  here  are  suggestive,  not  exhaustive.  Criteria  that  capture 
better  the  ideas  of  symmetry  and  size  might  be  developed  and  criteria  that  reflect 
other  goals  may  certainly  be  used. 


3.3  Random  Starts  for  the  Optimization  Procedure 

Multivariate  optimization  procedures  often  terminate  at  local  optimum,  rather 
than  at  the  global  optimum.  Because  of  this  phenomenon,  it  is  customary  to  start 
the  multivariate  optimization  procedure  with  many  different  initial  guesses  for  the 
solution  vector.  A  simple  method  of  obtaining  different  initial  values  for  the  vector 
0  is  to  generate  each  0^  as  a  uniform  random  variate  over  the  interval  (0, 2n).  This 
method  is  adequate  for  practical  purposes  but  it  is  subject  to  criticism  on 
theoretical  grounds.  The  columns  (or  rows)  of  G  can  be  interpreted  as  points  on 
the  surface  of  a  sphere  of  radius  1;  the  use  of  uniformly  distributed  0,;'s  does  not 
generate  a  sequence  of  random  orthogonal  matrices  whose  columns  are  uniformly 
distributed  over  a  sphere.  To  obtain  random  orthogonal  matrices  whose  columns 
represent  points  uniformly  distributed  over  a  sphere,  it  is  necessary  that  the  angles 
Oij  have  a  distribution  that  depends  on  the  value  of  j—i.  Anderson,  Olkin,  and 
Underhill  (1987)  give  the  following  method  of  obtaining  the  required  angles. 

1/2  y-*+i 

Let  $ij  =  arceos  {tyJi),  where  */_,  =  £  um  /  £  um 

m=l  m=l 

and  the  u’s  are  independent  random  chi-square  variates  with  one  degree  of  freedom. 
Anderson  et  al  (1987)  point  out  that  after  setting  tx  =  ux/(u  j+u2),  it  is  possible  to 
use  the  denominator  of  tx  as  the  numerator  in  f2  =  (u1+u2)  /  (u  i+u2+u3)  because 
the  properties  of  chi-square  variates  make  tx  and  t2  independent.  Thus,  for  fixed  t, 
the  k—i  Oij' 8  can  be  created  from  k—i+1  independent  chi-square  variates. 


4.  EXAMPLE 

Hardin  and  Sloane  (1992)  give  a  second- order  response  surface  design  for  four 
factors  in  16  runs  by  pairs  of  complex  numbers.  I  call  their  design  the  complex 
number  design.  Table  5  gives  the  real-valued  coordinates  for  the  four  factors  of  the 
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Table  5.  The  Complex  Number  Design 


Point 

** 

*4 

1 

0. 

0. 

0. 

0. 

2 

1. 

0. 

0. 

0. 

3 

-.5 

.866 

0. 

0. 

4 

-.5 

—.886 

0. 

0. 

5 

0. 

0. 

1. 

0. 

6 

0. 

0. 

-.5 

.866 

7 

0. 

0. 

-.5 

-.866 

8 

-.707 

0. 

-.707 

0. 

9 

-.707 

0. 

.354 

-.612 

10 

-.707 

0. 

.354 

.612 

11 

.354 

-.612 

-.707 

0. 

12 

.354 

-.612 

.354 

-.612 

13 

.354 

-.612 

.354 

.612 

14 

.354 

.612 

-.707 

0. 

15 

.354 

.612 

.354 

-.612 

16 

.354 

.612 

.354 

.612 

complex  number  design.  Rows  2-4  in  Table  5  correspond  to  the  pair  of  complex 
numbers  (wr,0)  for  r  =  0,  1,  and  2,  where  w=«xp(2jri/3);  rows  5-7  are  (0,u; r)  for 
r  =  0,  1,  2;  and  rows  8-16  are  2-1/2  (-wr,  -u*)  for  r  =  0, 1,  2  and  a  =  0,  1,  2  at  each 
value  of  r. 

The  design  has  been  oriented  to  minimise  the  number  of  levels  of  the  factors, 
but  Hardin  and  Sloane  (1992)  do  not  discuss  how  this  orientation  was  obtained. 
Notice  that  the  first  and  third  factors  cover  the  interval  —.707  to  1,  whereas  the 
second  and  fourth  factors  cover  the  interval  —.866  to  .866.  The  objectives  in 
rotating  this  design  were  to  obtain  (a)  a  symmetric  range  for  each  factor,  (b)  the 
same  range  for  each  factor,  (c)  the  same  set  of  levels  for  each  factor,  and  (d)  a  large 
region.  It  was  not  clear  if  all  these  objectives  could  be  attained,  nor  was  it  clear 
how  such  objectives  could  be  pursued  without  the  method  developed  in  this  report. 

The  first  step  in  exploring  this  problem  was  to  generate  100  random  0's  by 
the  method  of  Anderson  et  ai  (1987).  The  criteria  discussed  in  Section  3  were 
calculated  for  the  100  random  orientations  of  the  design  given  by  this  set  of  0's. 
Table  6  gives  the  mean,  standard  deviation,  minimum  and  maximum  of  the  criteria 
for  the  100  random  orientations,  for  the  design  as  initially  given  by  Hardin  and 
Sloane  (1992),  and  for  the  optimized  design.  For  the  optimization,  objective 
function  was  a  weighted  average  of  the  four  criteria;  the  weights  were  proportional 
to  the  reciprocals  of  the  standard  deviations  of  the  criteria.  Other  methods  of 
weighting  the  criteria  could  be  used.  For  example,  a  different  set  of  randomly 
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Table  6.  Criteria  for  Rot  at  ion*  of  the  Complex  Number  Design 


100  Random  Orientations 


Criterion 

Mean 

St.  Dev. 

Min. 

Max. 

Initial 

Optimised 

K|U 

1.78 

.04 

1.69 

1.84 

1.73 

1.71 

Aft 

.12 

.05 

.02 

.25 

.02 

.00 

ft#f» 

.40 

.13 

.07 

.67 

.59 

.00 

A  SAV 

.35 

.18 

.07 

.81 

.84 

.00 

varying  weights  could  be  used  with  each  starting  point  0. 


The  Nelder-Mead  downhill  simplex  method,  as  described  by  Press,  Flannery, 
Teukolsky,  and  Vetterling  (1986),  was  used  for  the  optimization.  Again,  this  was 
merely  a  convenient  choice;  other  optimization  methods,  such  as  direction-set 
methods,  annealing  methods  (see  Press  et  al  1986  for  a  description  of  these),  or 
genetic  algorithms  (see  Goldberg  1989)  could  be  used.  The  100  random  O' 8  were 
used  as  starting  points  for  the  optimization.  From  an  initial  point  0o,  the  downhill 
simplex  method  creates  a  simplex  by  generating  m~k(k—l)/2  additional  points  as 
follows.  For  i=l,  2, ...,  m,  let  0,  =  0o+^«i»  where  e,-  is  the  elementary  vector  with  a 
one  in  the  ith  place  and  zeros  elsewhere.  For  the  initial  simplex,  A  was  set  at  three 
radians.  Tne  downhill  simplex  method  moves  the  simplex  away  from  the  worst 
point  and  shrinks  the  simplex  until  the  relative  difference  between  the  objective 
function  at  the  best  point  and  the  objective  function  at  the  worst  point  reaches  a 
predetermined  value,  or  until  a  maximum  number  of  iterations  has  been  reached. 
The  termination  criterion  was  a  relative  tolerance  of  10-7  in  the  objective  function, 
or  a  limit  of  500  iterations.  Of  the  100  optimizations,  only  11  were  terminated  by 
the  500-iteration  limit.  The  angles  in  the  optimized  O' 8,  the  design  levels  in  W,  and 
the  four  criteria  for  a  design  were  generally  accurate  to  two  or  three  decimal 
places — as  judged  by  the  agreement  between  different  optimizations  that 
terminated  at  the  same,  or  equivalent,  rotations.  A  rerun  of  the  optimizations,  with 
the  optimized  O' 8  as  the  starting  points  and  A  =  .01,  added  another  decimal  place  of 
accuracy. 

Two  symmetric  orientations  were  found  that  have  the  same  value  by  each 
criterion.  The  two  orientations,  after  scaling  the  factor  levels  to  cover  the  interval 
—1  to  1,  are  given  in  Tables  7  and  8.  The  first  orientation  (Table  7)  can  be 
obtained  by  using  O'  =  (5x/6,  3.486,  5.161,  3tt/2,  2.299,  7tt/4),  changing  the  sign  of 
the  first  factor,  and  scaling  the  design  to  the  range  —1  to  1.  The  second  orientation 
(Table  8)  can  be  obtained  by  using  O'  —  (— x/6,  0,  arctan(21/2),  ?r/2,  3?r/4,  tt/4), 
changing  the  sign  of  the  fourth  factor,  and  scaling.  The  second  orientation  is 
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Table  7.  A  notation  of  the  Complex  Number  Design 


Point 

«» 

». 

*4 

1 

0 

0 

0 

0 

6 

8 

g 

0 

0 

4 

0 

0 

g 

g 

7 

-1 

b 

-d 

d 

5 

b 

-1 

d 

— d 

3 

d 

— d 

-1 

b 

2 

— d 

d 

b 

-1 

10 

1 

-d 

b 

d 

11 

-d 

1 

d 

b 

13 

d 

b 

1 

-d 

8 

b 

d 

-d 

1 

16 

f 

—a 

— C 

-h 

14 

—a 

f 

-h 

— c 

12 

-h 

— c 

f 

—a 

9 

— c 

-h 

—a 

f 

15 

—t 

— e 

—t 

— € 

NOTE:  a  =  3/2  1/5  -2«.121 

6  =  2^  a «. 172 
e  =  1-2-1/7  »  .293 
d  =  2*^-l«.414 
«  =  2-2‘/J  » .586 
/=2~‘/s  » .707 
j  =  2d  «  .828 
A  =  3c  w  .879 

preferable  to  the  first  because  it  has  an  additional  type  of  symmetry  not  found  in 
the  first  orientation:  if  a  factor  occurs  at  level  o,  it  also  occurs  at  level  —a. 


5.  ADDITIONAL  TOPICS 


5.1  Multiple  Solutions 

Any  angle  in  0  may  be  replaced  by  $ij  +  2n  without  changing  the  rotated 
design.  In  addition  to  the  multiple  solutions  created  by  the  periodic  nature  of  the 
sine  and  cosine  functions,  there  are  other  solutions  (0's)  that  represent  all 
permutations  of  the  columns  of  W  and  usually  some  permutations  of  the  rows  of  W. 
The  different  row  orders  that  can  be  obtained  by  rotation  are  called  equivalent  row 
orders;  the  implications  of  nonequivalent  row  orders  are  discussed  in  the  next 
paragraph. 
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Table  8.  Another  Rotation  of  the  Complex  Number  Design 


Point 

■i 

*» 

». 

•« 

1 

0 

0 

0 

0 

6 

g 

g 

0 

0 

4 

0 

0 

g 

g 

7 

-g 

0 

e 

— e 

5 

0 

-g 

— t 

e 

3 

— e 

e 

-g 

0 

2 

e 

— e 

0 

-g 

9 

-1 

-b 

0 

e 

12 

-b 

-1 

e 

0 

16 

e 

0 

-1 

-b 

14 

0 

e 

-b 

-1 

13 

1 

-d 

b 

d 

8 

-d 

1 

d 

b 

11 

d 

b 

1 

-d 

10 

b 

d 

-d 

1 

15 

—t 

— e 

— e 

— € 

NOTE:  b,  i,  e,  and  g  are  defined  in  Table  7. 


5.2  Proving  That  Two  Designs  Are  Rotations 

It  may  be  of  interest  to  show  that  two  Nxk  design  matrices  A  and  B  are,  or 
are  not,  rotations  of  each  other.  Applying  the  theory  of  generalized  inverses  to 
solve  A  =  BP  for  P,  with  the  idea  of  showing  that  P  is  an  orthogonal  matrix,  yields 
P  =  (B'B)_1B'A.  But  P  may  not  be  an  orthogonal  matrix  even  if  one  design  is  a 
rotation  (or  a  rotation  followed  by  a  reflection)  of  the  other.  Nonequivalent  row 
orders  for  A  and  B  will  make  P  a  nonorthogonal  matrix,  even  if  A  and  B  are 
rotations  of  each  other. 


5.3  The  Shape  of  the  Surface  Over  6 

If  one  considers  the  values  of  the  objective  function  as  a  surface  over  the 
possible  values  of  0,  then  the  shape  of  the  surface  depends  not  only  on  the  objective 
function  (the  criteria  used,  and  their  weighting),  but  also  on  the  orientation  of  the 
design  in  X.  Multivariate  optimization  procedures  are  sensitive  to  the  shape  of  the 
surface  over  the  parameter  space;  some  surfaces  make  the  multivariate  optimization 
problem  easy,  and  others  make  the  problem  difficult.  It  is  not  known  if  the 
optimization  process  can  be  improved  by  updating  the  design  in  X  whenever  a 
better  orientation  is  found. 
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5.4  Decomposition  of  P;  Range  of  0 

Anderson  et  al  (1987)  give  a  method  of  decomposing  an  orthogonal  matrix  P 
into  a  set  of  angles  6  that  can  be  used  to  construct  G  and  a  diagonal  matrix  D;  the 
diagonal  matrix  has  l’s  and  —  l’s  on  the  diagonal.  Let  P  =  GD;  then  G-1P  =  D. 
Because  G  is  the  product  of  k(k— 1)/2  matrices,  the  decomposition  of  P  into  D  and 
a  set  of  angles  is  performed  in  k(k—l)/2  stages;  each  stage  calculates  6  ij  and  then 
premultiplies  P  by  Gjy1.  The  0,y'«  are  extracted  in  the  standard  order  (1,2),  (1,3), 
...  (l,Jfc),  (2,3),  etc.  Let  d,y  =  aretan (p;,/p  „)  if  Pa> 0  and  let  0,y  =  tt/2  if  p  „  =  0. 
After  extracting  0^,  P  is  modified  by  using  the  equations 

P  *in  =  cos  (0,j)  pin  +  sin  {Oij)  Pjn  (10) 

and 

P*jn  =  -«»  (^»y)  Pin  +  COS  (Oij)  py„  (11) 

where  n=l,2, . . .  ,Jb  is  the  column  index  and  p*,n  and  p*yn  are  the  values  in  the  «th 
and  jth  rows  after  the  transformation.  Equations  (10)  and  (11)  are  similar  to 
equations  (4)  and  (5)  and  can  be  programmed  in  a  similar  manner.  At  the  end  of 
the  process  P  has  been  reduced  to  D. 

The  decomposition  of  P  given  by  Anderson  et  al  (1987)  always  yields  angles 
0,y  in  the  range  — tt/2  to  tt/2,  whereas  the  angles  in  the  multivariate  optimization 
procedure  are  allowed  to  take  on  any  value. 


5.5  Test  Cases  and  Results 

For  two  factors,  the  rotations  of  the  central  composite  design  and  of  the 
hexagon  design  given  in  Section  2  may  be  used  to  test  optimization  routines  and  the 
programming  of  equations  (4)  and  (5) .  For  k>  2  and  k  even,  the  central  composite 
designs  can  be  rotated  to  have  a  smaller  range  for  the  coded  factors  (and 
subsequently  scaled  to  cover  a  larger  region)  by  using  0,y  =  7r/4  for  <  odd  and 
j  = » +1,  and  Oij  =  0  otherwise. 

Doehlert’s  (1970)  uniform  shell  designs  for  k>  2  can  be  rotated  to  give  the 
factors  equal  ranges  by  using  0^  =  arcsin  [(j — *  +1) -1  ]  if  j<k  and 
6^  =  —  arenn\{j— i+l)-1/2  ]  if  j  =  k.  The  designs  produced  by  this  rotation  have 
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seven  levels  of  the  factors,  except  for  the  three-factor  design  (three  levels — the 
Box-Behnken  design)  and  the  eight-factor  design  (live  levels).  It  is  necessary  to 
change  the  sign  of  the  last  factor,  as  the  preferred  form  for  the  new  orientation  is 
both  a  rotation  and  a  reflection  of  the  uniform  shell  designs.  The  designs  for  even  k 
can  be  orthogonally  blocked  in  the  same  manner  as  the  design  in  Table  4. 


6.  SUMMARY 

Figures  illustrating  rotations  of  the  hexagon  design  and  the  octagon  (or 
central  composite)  design  showed  how  rotations  can  produce  a  symmetric  treatment 
of  the  factors,  a  larger  experimental  region,  or  both.  A  rotation  of  the  eight-factor 
uniform  shell  design  (with  levels  —1,  —.5,  0,  .5,  and  1  for  each  factor)  was  given. 

Rotations  of  a  response  surface  design  can  be  obtained  from  standard 
multivariate  optimization  techniques.  This  method  requires  parameterization  of 
orthogonal  matrices  by  a  set  of  parameters  that  can  be  varied  independently  and 
the  specification  of  the  objective  function  for  the  multivariate  optimization  routine. 
The  factorization  of  an  orthogonal  matrix  into  matrices  that  represent  planar 
rotations  was  used  to  represent  the  orthogonal  matrix  by  a  set  of  planar  angles  that 
can  be  varied  independently.  For  the  objective  function,  design  criteria  that 
represent  the  size  of  the  design  region  and  the  symmetry  of  the  design  matrix  were 
given.  The  specific  criteria  developed  were  the  maximum  range  of  the  factors,  the 
equality  of  the  factor  ranges,  the  symmetry  of  the  factor  ranges  about  the  center 
point  of  the  design,  and  the  equality  of  the  sets  of  factor  levels. 

Multivariate  optimization  procedures  require  a  starting  value  (initial  guess) 
for  the  vector  of  planar  angles  that  represent  an  orthogonal  matrix.  The  method  of 
Anderson  et  al  (1987)  for  generating  the  angles  will  produce  a  sequence  of  random 
orthogonal  matrices  that  represent  rotations  uniformly  distributed  over  a  sphere. 

The  method  for  rotating  a  design  was  applied  to  the  four -feu: tor  complex 
number  design  of  Hardin  and  Sloane  (1992).  The  criteria  for  the  desirability  of  a 
rotation  were  calculated  for  100  random  orientations;  the  criteria  were  then  divided 
by  their  standard  deviations  to  form  a  weighted  average  for  use  as  the  objective 
function  for  the  multivariate  optimization.  Two  symmetric  orientations  of  the 
design  were  found. 
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Appendix:  Proof  of  Invariance  of  D-  and  G- Efficiency  Under  Rotation 

The  proofs  for  a  second-order  response  surface  design  are  more  difficult  than 
the  proofs  for  a  first-order  model.  I  will  give  only  the  proofs  for  the  second-order 
case  in  detail;  the  proofs  for  the  first-order  case  are  similar  to  those  for  the  second- 
order  case.  These  proofs  show  that  the  D-  and  G -efficiencies  of  a  second-order 
response  surface  design  are  not  changed  by  a  rigid  rotation  of  the  coordinate  axes. 
In  a  spherical  region  rotation  of  the  coordinate  axes  is  equivalent  to  a  rotation  of 
the  design  points.  For  a  cuboidal  region,  rotation  of  the  coordinate  axes  is 
equivalent  to  rotation  of  the  design  points  and  the  cuboidal  region  simultaneously 
(so  that  the  design  points  maintain  the  same  relationship  to  the  region  boundaries 
as  they  had  before  the  rotation). 

Write  the  first-  or  second-order  model  as 


y  =  M,  /3  +  e  (A.l) 

where  y  is  the  iVxl  vector  of  responses,  M,  is  the  Nxp  model  matrix,  /?  is  the  pxl 
vector  of  coefficients,  and  e  is  the  Nx  1  vector  of  errors.  A  rotation  (possibly 
followed  by  a  reflection)  of  the  Nx  k  design  matrix  X  is  given  by 

W  =  XP,  (A.2) 

where  P  is  an  orthogonal  matrix.  The  model  matrix  for  the  rotated  design  is 
and  the  coefficient  vector  for  the  rotated  design  is  denoted  7. 


Let  Mw  =M,  Q.  For  a  first-order  model, 


Q  = 


1  o' 

0  p 


(A.3) 


Q  is  therefore  an  orthogonal  matrix.  For  the  second-order  model  the  matrix  that 
maps  Mz  to  Mw  is  not  an  orthogonal  matrix. 

Reparameterize  the  model  (1)  for  the  second-order  case  by  replacing  the 
interaction  (product)  terms  z,  z;-  by  21/2  z,-  z;-  and  the  interaction  parameters  /?,;  by 
2-1/2  0ij.  To  write  the  reparameterized  model  in  matrix  notation,  I  introduce  a 
diagonal  matrix  C  that  has  l’s  in  the  positions  corresponding  to  the  constant  term, 
linear  terms,  and  squared  terms  and  21/2  in  the  positions  corresponding  to  the 
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interaction  terms.  The  reparameterized  model  is  written 
y  =  M,CC~1/3  +  c. 


(A. 4) 


After  rotation  of  the  design  matrix  X  by  P,  the  model  (4)  becomes 

y  =  MwCC“S+c-  (A.5) 

For  the  reparameterized  models, 

M*C  =  MJCQ,  (A.6) 

1  O'  0' 

where  Q  =  0  P  0 

0  0  pl2! 

pt2l  is  the  second  Schlaflian  matrix  derived  from  P  and  maps  the  second-order 
terms  of  the  reparameterized  model  to  the  second- order  terms  of  the  rotated, 
reparameterized  model  (see  Box  and  Hunter  1957).  P^  is  an  orthogonal  matrix 
and  therefore  Q  is  also  an  orthogonal  matrix. 

To  show  the  invariance  of  D -efficiency,  I  need  to  show  that  the  determinant 
of  M^MW  is  equal  to  the  determinant  of  MjM,.  Using  the  elementary  properties 
of  matrices  and  determinants, 

det [(M„C) '  MWC]  =  dct(C'M;MwC)  =  det  (C')  det(M;,M„)  det( C).  (A.7) 

Applying  equation  (A.6)  to  the  left  hand  side  of  equation  (A.7)  gives 
det[(MtBC)'MM,C]  =  det  [(M^QJ'MjCQ]  =  dct  (Q'C'M;MxCQ)  = 
det(Q')det(C')det(M;M,)def(C)det(Q)  =  det  (C')  det  (M;Ms)  det  (C)  (A.8) 

because  det(Q)  =  ±1.  Combining  equations  (A.7)  and  (A.8)  gives 

det  (C')  det  (M;  M„)  det  (C)  =  det  (C')  det  (Mi  M,)  det  (C).  (A.9) 

Because  det  (C)  **=(), 
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(A.10) 


To  show  the  invariance  of  G-efficiency  under  rotation,  I  need  to  show  that 
the  maximum  prediction  variance  within  the  experimental  region  remains 
unchanged  when  the  design  is  rotated.  To  do  this,  I  will  show  that  the  prediction 
variance  at  every  point  in  the  region  remains  unchanged  when  the  design  is  rotated. 
Let  x'  and  w'  be  lxp  vectors  that  represent  the  same  point  in  the  models  Mz  and 
Mu  respectively.  In  the  reparametrized  models  using  MZC  and  Mu  C,  the  point  is 
represented  by  w'C  and  x'C,  respectively.  Therefore  w'C  =  x'CQ  and  (its 
transpose)  C'w  =  Q'C'x. 

Now  examine  the  effect  of  the  reparameterization  on  the  prediction  variance: 

For(y \w‘C)/a2  =  w'C[(M„C)'  M.Cj^C'w  =  w'C[C'M;MliC]-1C'w  = 

w'CC~1(M^Mt0)_1C'_1  C'w  =  w'(M^Mtt)-1w  =  Var(y  |  w')/o  2.  (A.ll) 

Next,  examine  the  effect  of  the  rotation  on  the  prediction  variance  of  the 
reparameterized  model: 

Var{y  |  w'C)/o  2  =  w'C[(M„C)' M^Cj^C'w  = 
x'CQ[(MzCQ)'(MzCQ)]-1Q'C'x  =  x'CQlQ'C'MiM^CQl^Q'C'x  = 
x'cqq-1c-1(m;mx)-1c'-1q'-1qc'x  = 

x'(M;Mz)"1X=  Vdr(y  |x')/<7  2.  (A.12) 

Combining  equations  (A.ll)  and  (A.12)  yields  Var{y  Jw')  =  For  (y  |x'). 
APPENDIX 
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